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The chemical disorder intrinsic to high entropy alloys inevitably creates diffuse scattering in 
their x-ray or neutron diffraction patterns. Through first principles hybrid Monte Carlo/molecular 
dynamics simulations of two BCC high entropy alloy forming compounds, CrMoNbV and NbTiVZr, 
we identify the contributions of chemical disorder, atomic size and thermal fluctuations to the diffuse 
scattering. As a side benefit, we evaluate the reduction in entropy due to pair correlations within 
the framework of the cluster variation method. Finally, we note that the preference of Ti and Zr 
for HCP structures at low temperature leads to a mechanical instability reducing the local BCC 
character of NbTiVZr, while preserving global BCC symmetry. 

PACS numbers: 


I. INTRODUCTION 

While the Bragg peaks of diffraction patterns reveal the global symmetries and average structures of crystals, diffuse 
scattering intensity located between the peaks contains information about the deviations from the average. Many 
factors create diffuse scattering: chemical disorder creates smooth diffuse background due to contrasting scattering 
form factors; variation in atomic sizes and thermal motion (phonons) leads to diffuse wings that diverge in intensity 
near the Bragg peaks each with its own characteristic angular variation; a variety of extrinsic defects each with their 
own diffuse signature. 

Because each type of disorder creates diffuse intensity of specific and quantifiable intensity patterns, a great deal 
of physical information can be extracted from diffuse scattering experiments. Disorder in atomic form factors creates 
a uniform diffuse background, but short-range chemical order can cause the background to concentrate in specific 
regions between or away from Bragg peaks depending on whether the interatomic interactions favor local association 
of like or unlike chemical species. Thermal diffuse scattering is dominated by long wavelength phonons and thus the 
bulk, shear and other elastic moduli can be extracted from its angular variation. Atomic size scattering arises from 
long-range strain fields caused by mismatch of an atoms size to its environment. Its strength depends on the mismatch 
and on particular combinations of elastic constants. 
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We explore the diffuse scattering patterns of two high entropy alloys by means of computer simulation. The 
compounds NbTiVZr l| and CrMoNbV 2j are drawn from refractory metals of, respectively, the fourth and fifth, 
and the fifth and sixth columns of the periodic table. Thus the atomic form factors and the contrasts among the form 
factors are similar for the two cases. Because both belong to squares of the periodic table 2, 3}, we anticipate the 
greatest contrast in size to lie along the positive diagonals of the squares, namely Zr-V and Nb-Cr. One interesting 
point about these alloy systems is that all constituent elements of CrMoNbV individually form body centered cubic 
(BCC) lattices, while two of the elements (Zr and Ti) of NbTiVZr are BCC only at high temperature but prefer 
hexagonal (HEX) lattices at low temperature. Indeed, the BCC structures of these elements are mechanically unstable 
at low temperature. 

In order to efficiently explore the configurational ensemble in the solid state, we utilize a hybrid MC/MD method 
that alternates Monte Carlo swapping of atomic species with conventional molecular dynamics 3|. A side benefit of 
this simulation method is that we obtain partial pair correlation functions of the alloy systems. From these correlation 
functions we can see the impact of atomic size differentials and preferences in short-range chemical order. Remarkably, 
we see a loss of short-range BCC geometry in the case of NbTiVZr. This results in stronger diffuse scattering in 
NbTiVZr than in CrMoNbV. The existence of short-range chemical order requires that the actual configurational 
entropy lie below its ideal value of fee log (4). We apply the formalism of the cluster variation method (CVM) to 
estimate the actual entropy (at T=1200K) of NbTiVZr as kg log (3.93) and CrMoNbV as kg log (3.76). 


II. RESULTS 


A. Swap rates and pair correlations 


Using hybrid MC/MD we pre-anneal structures of N = 128 atoms (32 of each species) at T=1800K and then 1500K 
prior to lengthy data collection runs at 1200K. Our data collection runs extend for 15ps of molecular dynamics and 
1500 attempted swaps for NbTiVZr, and 28ps of MD and 2800 attempted swaps for CrMoNbV. Acceptance rates for 
swaps indicate the similarity of chemical species in size and bonding preferences. Characteristically, the swap rates 
are lowest across the positive diagonal of periodic table squares 2] (i.e. Nb with Cr, and Zr with V) and highest along 
the negative diagonal (i.e. V with Mo, and ti with Nb). The relatively large acceptance rates for most pairs indicate 
the efficiency of our MC/MD method and suggest a likely large entropy of mixing. 


Pair correlation functions are more revealing at lower temperature where atomic motion is reduced. To mimic an 
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Swap 

NN Bonds 

NNN Bonds 

a\/3 

Zr Ti Nb V 

Zr Ti Nb V 

Zr Ti Nb V 

Zr 

- 0.19 0.24 0.06 

0.055 0.067 0.060 0.068 

0.073 0.060 0.064 0.053 

Ti 

- 0.57 0.19 

0.058 0.065 0.061 

0.061 0.064 0.066 

Nb 

- 0.15 

0.061 0.064 

0.055 0.067 

V 

- 

0.057 

0.064 


TABLE I: Monte Carlo swap rates and bond counts for Nb-Ti-V-Zr quaternary at T = 1200K. Bond counts N a> p count number 
of /3-type neighbors of atom type a , where a labels rows and /3 labels columns. Elements are arranged in order of decreasing 
BCC lattice constant. 



Swap 

NN Bonds 

NNN Bonds 

q\/3 

Nb Mo V Cr 

Nb Mo V Cr 

Nb Mo V Cr 

Nb 

- 0.44 0.23 0.05 

0.049 0.063 0.062 0.076 

0.080 0.067 0.056 0.045 

Mo 

- 0.50 0.21 

0.058 0.067 0.062 

0.061 0.064 0.059 

V 

- 0.45 

0.055 0.065 

0.066 0.066 

Cr 

- 

0.048 

0.080 


TABLE II: Monte Carlo swap rates and bond statistics for Cr-Mo-Nb-V quaternary at T = 1200K. Other details as in Table [I] 


experiment where the alloy is formed at high temperature then cooled down, we quench from T=1200K down to 300K 
by performing molecular dynamics alone without Monte Carlo swapping. The chemical order remains characteristic of 
T=1200K even while the phonons equilibrate at 300K, as is the case in typical metallurgical experiments. Inspecting 
Fig. Q] note that CrMoNbV exhibits more ideal BCC structure at short range, with the nearest neighbor (NN) and 
next-nearest neighbor (NNN) peaks clearly resolved, than is the case for NbTiVZr, where these peaks are poorly 
resolved for Nb and V, and not resolved at all for Zr and Ti. Nonetheless, the long-range BCC structure is well 
preserved in both compounds as can be seen in the further neighbor peaks. A similar local instability was observed 
in ternary HfNbZr ^ which likewise mixes elements from the fourth and fifth periodic table columns (be. HEX and 
BCC ground states). 

Because BCC is a loosely-packed structure, instantaneous atomic positions can be mapped back to the ideal BCC 


lattice sites that they are closest to. This allows us to uniquely identify NN and NNN bonds and report their statistics, 
as is done in Tables U and |TI1 Notice that the largest number of NN pairs occurs between the elements on the positive 
diagonals of the squares (Zr-V and Nb-Cr). This is reflected in Fig. [Tj in the respective strengths of the peaks. 
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FIG. 1: Pair distribution functions of (a) NbTiVZr and (b) CrMoNbV at T=300K quenched from 1200K. Each panel shows 
the four partial pair correlation functions for the element named above. The partials are color coded, e.g. in (a) for NbTiVZr, 
black indicates Ti, red indicates V, green indicates Zr and blue indicates Nb. Thus under Ti the black curve is Ti-Ti and the 
red curve is Ti-V, etc.. Bars at top indicate the corresponding correlations in the pure BCC element, e.g. elemental Ti has 8 
neighbors at r = 2.9 A, 6 at 3.3, 12 at 4.7, 24 at 5.5 and 8 at 5.7. 


B. Diffuse scattering 


An instantaneous configuration is described by positions r m that are displaced from the ideal lattice sites R m . We 
specify chemical occupation at position r m by c“ defined as 1 if occupied by species a and 0 otherwise. The scattering 
amplitude at reciprocal space position Q is [5-8] 


^(Q) = EE/«i e " Q ' rm ' (!) 

m ol 

and its square is the scattering intensity /( Q) = |A(Q)| 2 . Here f a is the scattering form factor f a . For the present 
analysis we utilize x-ray form factors which we approximate as the atomic number, f m ss Z mi independent of Q. Note 
that the set of atomic numbers present in NbTiVZr (22, 23, 40 and 41) exhibits contrast similar to the set present 
in CrMoNbV (23, 24, 41, and 42). The natural logarithm of the intensity is plotted in Fig. (2j which is measured in 
reciprocal lattice units (RLU), (27r/a), where a is the conventional cubic lattice constant. Bragg peaks of the BCC 
lattices occur at Q = G = (27r/a)(Ff, K, L) where H 1 K and L are integers such that H + K + L is even. For Q close 
to a Bragg peak G we shall denote the deviation by k = Q — G. 

The average scattering intensity /( Q) = (|A(Q)| 2 ) separates into the Bragg component of the average lattice 
Ib( Q) = |(A(Q))| 2 , and the diffuse part due to fluctuations Id(Q) = (|^4(Q) — (A(Q))| 2 ) = I — I B . Owing to 
configurational fluctuations we take averages over both atomic displacements u m = r m — R m and chemical species 
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fluctuations Ac“ = c“ — c“, yielding [7| 


/d(Q) = ^ ACAc£e- ,< *< u ’»»>) 


( 2 ) 


a/3 


where we introduce the relative coordinates R mn = R n — R m and u mn = u„ u m . For sufficiently small Q • u, we 
may expand the exponential 


- iQ - (Umn ) wl-iQ- (u m „) - i(Q • (u mn )) 2 + • • • 


( 3 ) 


Separating the diffuse scattering into components, we set Id = To + Ji + Ii + • • •, where each term arises from 
successively higher powers of Q ■ u. Here Iq reflects short-range chemical order, 7i is known as the size effect, and I 2 
includes both the thermal diffuse scattering due to lattice vibrations (phonons) and also the Huang scattering due to 
correlated lattice strains caused by differing atomic sizes. 

To model chemical disorder, assume that all atoms sit on BCC lattice sites {Rm} with chemical species randomly 
distributed. For Q not on a Bragg peak, i.e. k = Q - G ^ 0, the sum over m and n in Eq. d2J) simplifies to 
Iq/N ~ (/ 2 ) — (/) 2 independent of k, resulting in a uniform diffuse background. However, our neighbor statistics given 
in Tables Q] and El show the distributions of species is not fully random, but rather exhibits short-range correlations. 
Unlike species (e.g. Zr with V in NbTiVZr, or Nb with Cr in CrMoNbV) prefer to occupy nearest neighbor bonds. 
On a BCC lattice such sort-range chemical order can lead to partial B2-like ordering in which broad diffuse maxima 
arise at previously forbidden positions where H + K + L is odd, as indeed we see in the case of NbTiVZr (Fig. [5}a)). 
Despite the short-range order present in CrMoNbV, we see no B2-like peaks in its diffraction pattern, which could 
indicate more complex order such as the B 23 predicted in NbMo [9|. 

In fact, the atoms in our models do not sit at ideal lattice sites. The random distribution of unlike atomic sizes 
cause substantial relaxations off the lattice sites even at T=0K, and in addition they are subject to ordinary thermal 
fluctuations. Atomic size effects enter both at first and second order. At the first order the “atomic size effect” 


yields [8] 


h{ Q) = - 5Z faf(3 ' Umra ) sinQ ■ R„ 


( 4 ) 


a.f3 


Nonvanishing I\ depends on a correlation between the sign of u mn and the strengths of the form factors at sites m 
and n. The long-range displacement u rnn ~ ea/i ? 2 has Fourier transform u(k) ~ e^/k. Thus close to a Bragg peak 
at Q = G + k, the size effect scattering diverges [7} as 7i(Q) ~ — G • e^/k and will be asymmetric on the two sides of 
G. For example, if atoms with large / a tend to move apart and small f a tend to move closer, then I\ is positive for 


negative k and negative for positive k 8 |. This asymmetry can give the impression of a shift in the peak position 


id- 
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FIG. 2: Diffuse scattering patterns (natural logarithm of intensity J(Q)) of NbTiVZr (a, c and e) and CrMoNbV (b, d and f) 
at temperatures T = 300K (a and b) and T = OK (c and d). Patterns were obtained from samples of 256 atoms in 8 x 8 x 2 
supercells of the BCC unit cell. Curves marked “Lattice” in (H, 0, 0) scans (e and f) arise from placing the atoms at their ideal 
lattice positions. 

At second order ordinary thermal diffuse scattering due to phonons, and also Huang scattering due to atomic 
sizes, both contribute diffuse scattering that diverges as 1/k 2 in the vicinity of the Bragg peaks, each with its own 
characteristic angular distribution. Ordinary thermal fluctuations are characterized by phonons whose dispersion 
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relations are of the form u)j = Cjk , where j = 1..3 enumerate the independent modes of vibration whose polarization 


vectors are ej( k). The resulting diffuse intensity is proportional to 


ul 


(iG-u k ivE iG, tr ,r 


(5) 


As transverse sound speeds are typically less than longitudinal, the diffuse scattering will usually be stronger in 
directions k perpendicular to G than parallel to it. Such elongations can be seen in the diffuse pattern for CrMoNbV 
in surrounding the (200) and (400) peaks. Also visible are streaks running in [110] directions similar to those observed 


in other BCC metals 


12], These could indicate a low frequency phonon mode with a polarization component in the 


[110] direction. Notice that the streaks are stronger on the low Q side of the Bragg peaks (i.e. negative G • ek) than 
on the high Q side, presumably as a result of the asymmetry of the size effect scattering I\. In fact, there are local 
maxima of scattering in CrMoNbV that shift from, e.g. H = 3.1 at T = 300K to H = 3.5 at H = OK, as can be seen 
in Fig. [2f indicating the diminishing contribution of the thermal component of I 2 relative to I\. 

Huang scattering arises from long-range strains arising from atomic size mismatch and hence is temperature inde¬ 
pendent. Originally derived for dilute impurities, a similar effect is present in concentrated alloys such as high entropy 
alloys. Because these strains are primarily longitudinal, only the component of k parallel to G enters, and the diffuse 
intensity is proportional to 


which is equivalent to Eq. [5] setting ej parallel to k. The Huang scattering consists of figure-eights that are charac¬ 
teristically elongated in the longitudinal (radial) direction in reciprocal space. The low T patterns in Fig. [2] indeed 
show such elongation, and again the asymmetry in the sign of k arises from the contribution of I\. The figure-eight 
pattern is not clearly resolved in the case of NbTiVZr, due to transverse strain that might indicate the presence of 
multisite correlations. 


C. Configurational entropy 


We evaluate the configurational entropy using the approach of the cluster variation method (CVM ljj 14]) by 
starting with the mean field (Bragg-Williams) entropy, then modifying it through the inclusion of factors that reduce 
the entropy by correcting for local correlations. First we define numerical values for symbols associated with the 
empty lattice, isolated points (P), near-neighbor bonds (NN), next-nearest-neighbor bonds (NNN), triangles (TRI) 










NbTiVZr 

CrMoNbV 


4 3.971 3.931 (3.789) 
4 3.886 3.759 (3.657) 


method 


TABLE III: Entropies of NbTiVZr and CrMoNbV at T = 1200K in the mean field approximation (MF 151), the quasichemical 
approximation (QC Q ), an improved QC incorporating both NN and NNN correlations, and the cluster variation (CVM) 


171 ]. CVM entropies are poorly converged owing to limited sampling statistics of the four-point correlation. Values 


quoted are phase space volumes exp ( S/ks ) where S is the entropy per site. 


and tetrahedra (TET) as 


{1} = M\, = = (2MM)!, 


o(3 


{—} = (w a pM)\, {.A,} = Y[{v a ^M)\, {<£>} = ( z a 0 jS M )\, 

a/3 o(3~y a^'yS 


(7) 


( 8 ) 


where M is the total number of lattice sites, x a are the fractions of species a , y a p is the frequency of NN bonds 
between species a and /3, etc. Next, we introduce combinatorial factors 


f!p = 


ill 

M’ 


&NN — 


M 8 




^NNN = 


M 6 


W 3 {1} 3 ’ 


Qtet = 


Note that jglogQp ss — J2 a l°g x a is the mean field entropy 




(9) 


15j. Multiplying flp by Qnn reduces the entropy 


by an amount related to the deviation of NN bond frequencies y a p from the uncorrelated frequency x a xp. This level 


of approximation also goes by the name “quasichemical approximation” 


16]. Multiplying instead by Qjvatjv would 


do the same for the NNN bonds, so the product Qnn^nnn given in Table Hill incorporates both effects, that is, it 
corrects the quasichemical approximation through inclusion of NNN correlations. Finally, including the remaining 
factor Qtbt results in the conventional CVM entropy expression for BCC lattices 17]. 


By accumulating cluster occupation statistics during our MC/MD simulation, we are able to evaluate the numerical 
factors. Table un presents the resulting entropies (in the form of the phase space volume per site, which is the 
exponential of the entropy per site) for our two compounds. Note that these are configurational entropies only [i.e. 
no vibrational or electronic contributions), and that they reflect the chemical disorder at T=1200K. The entropies 
decrease monotonically as additional correlations are included. However, the CVM entropy must not be taken as 
an upper bound, because our sampling statistics are insufficient for the four point correlation function Zap^g, and 
additional run time would increase this value. The proper values lie between NNN and CVM, and we recommend 
taking the NNN value as our best estimate. Interestingly, the entropy of NbTiVZr lies closer to the ideal (mean 
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field) value than CrMoNbV does. This might reflect the extreme short-range disorder of NbTiVZr, which seemingly 
precludes establishing strong local chemical order. 


III. CONCLUSION 

In summary, we have modeled the structures of two body-centered cubic refractory high entropy alloys, NbTiVZr and 
CrMoNbV. We discover a local instability of NbTiVZr that reduces the short-range BCC structure while preserving the 
long-range lattice. This instability is likely due to the preference of Ti and Zr for HEX structures at low temperature. 
Chemical disorder, atomic size differences and thermal fluctuations lead to diffuse x-ray scattering patterns, with 
chemical and size effects dominating at room temperature in the case of NbTiVZr, while thermal fluctuations dominate 
in the case of CrMoNbV. A hint of B2-like chemical order is evident in NbTiVZr but is surprisingly absent in 
CrMoNbV. A side benefit of our simulations is the ability to evaluate configurational entropy from the pair correlation 
functions. We find the entropies at T = 1200K are ks In (3.93) for NbTiVZr and ks In (3.76) for CrMoNbV. 
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